//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// g++ -I /home/tiago/astro_local/lib/ -O2 cur_luz.cpp -o cur_luz
//
// cur_luz.cpp - Nova versao do programa curluz2.c agora escrito em C++.
// Inclui o biblioteca de leitura de arquivos dataFileReader.h (Zabot A.)
// de modo que comentarios podem ser inseridos nos arquivos de entrada. Os
// comentarios devem ser precedidos de #.
// 
// Modo de utilizacao do programa:
// cur_luz arq_in arq_out
// onde arq_in o arquivo com os parametros de entrada e arq_out o arquivo 
// para serem escritos os resultados.
//
// Formato do arquivo de entrada:
// 1a. linha de dados (sem contar linha de comentario):
//   i  - inclinacao (graus).
//   q  - razao de massa (M2/M1).
//   f2 - fluxo da secundaria em fase 0. 
//   fd - fluxo do disco (constante somada a curva de luz).
//   S/N- Relacao Sinal-Ruido na curva
// 2a. linha de dados:
//   limb_law  - lei de obscurecimento de borda.
//   grav_coef - coeficiente do obscurecimento gravitacional.
//   Tpole     - Temperatura do polo da estrela.
//   nspot     - Numero de spots 
// 3a. linha de dados:
//   (coeficientes de obscurecimento de borda)
// 4a. a (4+nspot) linhas de dados: 
//   theta - posicao theta do spot
//   phi   - posicao phi do spot
//   dt    - largura a meia altura em theta
//   dp    - largura a meia altura em phi
//   Is    - Intensidade em relacao a I_L1
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include <iostream>
#include <cmath>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <cstring>
#include "light_curve.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"
//#include "light_curve.cpp"
//#include "roche.cpp"
//#include "lib_zabot/dataFileReader.h"

using namespace std;

int definicoes_obscurecimento(const int law);
void show_par(ostream& out, vector<vector<double> > &par);
double S(vector<double> &);
double S_gauss(vector<double> &,vector<double> &);

int main(int argc, char *argv[])
{
	// Testar argumentos de entrada...
	
	if(argc < 3){
		cout << "Modo de utilizacao do programa:\n\n"
		<< "\timOGL arq_p arq_out arq_img\n\n"
		<< "Programa que passa imagens .dat para .ogl\n"
		<< "Onde arq_p o arquivo com os parametros de entrada e arq_out o\n" 
		<< "arquivo para ser escrita a imagem final. Opcionalmente pode\n"
		<< "se incluir um arquivo com as fases desejadas para ser\n"
		<< "calculados o modelo ([fase.dat]).\n\n"
		<< "Formato do arquivo de entrada:\n"
		<< "1a. linha de dados (sem contar linha de comentario):\n"
		<< "  i  - inclinacao (graus).\n"
		<< "  M1 - Massa da Primaria.\n"
		<< "  M2 - Massa da Secundaria. (q = M2/M1)\n"
		<< "  P  - Periodo Orbital (hrs).\n"
		<< "2a. linha de dados:\n"
		<< "  f2 - fluxo da secundaria em fase 0.\n"
		<< "  grav_coef - coeficiente do obscurecimento gravitacional.\n"
		<< "  limb_law  - lei de obscurecimento de borda.\n"
		<< "  nspot     - Numero de spots.\n"
		<< "  Tpole     - Temperatura do polo da estrela.\n"
		<< "  fd - fluxo do disco (constante somada a curva de luz).\n"
		<< "3a. linha de dados:\n"
		<< "  (coeficientes de obscurecimento de borda)\n"
		<< "4a. a (4+nspot) linhas de dados:\n"
		<< "  theta - posicao theta do spot\n"
		<< "  phi   - posicao phi do spot\n"
		<< "  dt    - largura a meia altura em theta\n"
		<< "  dp    - largura a meia altura em phi\n"
		<< "  Is    - Intensidade em relacao a I_L1\n\n"
		<< "Tiago Ribeiro - OpenGL/LinuxC++ - 2007\n";
		exit(1);
	}
	
	// Input File Definition - IFD:  
	//     Number of fixed Parameters lines      -  IFD_NP
	//     Number of columns Parameters          - *IFD_P (array[IFD_NP])
	//     Spot Line(onde com. def. de spot - 1) -  IFD_SL
	//     Number of Spot Parameters             -  IFD_SP
	//     posicao da lei de obs. de Borda       -  IFD_B
	
	const int IFD_NP = 2;
	const int IFD_P[IFD_NP] = {4, 7};
	const int IFD_SL = 3;
	const int IFD_SP = 5;
	const int IFD_B = 2;
	
	DataFileReader fin(argv[1]);
	
	vector<vector<double> > in_par(IFD_NP+1);
	
	//Lendo parametros orbitais
	
	for(int i = 0; i < IFD_NP; i++){
		++fin;
		if(fin.numberOfColumns() == IFD_P[i]){
			double data_in;
			for(int j = 0; j < IFD_P[i]; j++){
				data_in = fin.readColumn<double>(j+1);
				if(i == 0 && j == 0) data_in *= M_PI/180;
				in_par[i].push_back(data_in);
			}
		}
		else { 
			ostringstream erro;
			erro << "Numero de parametros errado\n\tArquivo: " 
			<< argv[1] << "\n\tLinha  : " << fin.lineNumber();
			aborta(erro.str());
		}
		/*if(fin.eof()) {
		 ostringstream erro;
		 erro << "Numero de parametros errado\n\tArquivo: " 
		 << argv[1] << "\n\tLinha  : " << fin.lineNumber();
		 aborta(erro.str());
		 }*/
	}
	
	//Lendo coeficientes do obscurecimento de borda
	++fin;
	if(fin.numberOfColumns()==definicoes_obscurecimento(int(in_par[1][IFD_B]))){
		double data_in;
		for(int i = 0; i < definicoes_obscurecimento(int(in_par[1][IFD_B])); i++){
			data_in = fin.readColumn<double>(i+1);
			in_par[2].push_back(data_in);
		}
	}
	else {
		ostringstream erro;
		erro << "Numero de parametros errado\n\tArquivo: " 
		<< argv[1] << "\n\tLinha  : " << fin.lineNumber();
		aborta(erro.str());
	}
	
	
	/* if(fin.eof()) {
	 ostringstream erro;
	 erro << "Numero de parametros errado\n\tArquivo: " 
	 << argv[1] << "\n\tLinha  : " << fin.lineNumber();
	 aborta(erro.str());
	 }
	 */
	//Lendo parametros do spot
	
	for(int i = 0; i < int(in_par[1][3]); i++){
		++fin;
		if(fin.numberOfColumns() == IFD_SP ){
			double data_in;
			in_par.push_back(vector<double> ());
			for(int j = 0; j < IFD_SP; j++){
				data_in = fin.readColumn<double>(j+1) ;
				in_par[IFD_SL+i].push_back(data_in);
			}
		}
		else {
			ostringstream erro;
			erro << "Numero de parametros errado\n\tArquivo: " 
			<< argv[1] << "\n\tLinha  : " << fin.lineNumber();
			aborta(erro.str());
		}
		
	}
	
	fin.close();
	
	show_par(cout, in_par);
	
	surface SecIm;
	
	
	Roche Star2(4,in_par);
	//if(argc == 3) Star2.faz_super();
	//Star2.spot_maker_test();
	Star2.spot_maker();	
	surface Iref = Star2.Star(); 
	
	cout << "# - Model map entropy is, S = " << S(Iref.I) << endl;
	
	ofstream saida(argv[2]);
	
	if(!saida){
		ostringstream erro;
		erro << "Arquivo " << argv[2] << " nao pode ser aberto!";
		aborta(erro.str());
	}
	
	show_par(saida, in_par);
	
	for(int i = 0; i < Iref.theta.size(); i++){
				
		saida << Iref.theta[i] << " "
			<< Iref.phi[i] << " "
		<< Iref.px[i] << " "
		<< Iref.py[i] << " "
		<< Iref.pz[i] << " "
		<< Iref.nx[i] << " "
		<< Iref.ny[i] << " "
		<< Iref.nz[i] << " "
		<< Iref.I[i] << " "
		<< Iref.vx[i] << " "
		<< Iref.vy[i] << " "
		<< Iref.vz[i] << endl;

	}
	
	saida.close();
	
	cout << "--\n#-> Show rimg figures on python with:\n#\n"
		<< "import healpy as H\n"
		<< "import numpy as np\n" 
		<< "map = np.loadtxt('"<<argv[2]<<"').T[8]\n"
		<< "H.mollview(map,rot=(-90,-90),cmap=cm.hot)\n--\n\n";
	
	return 0;
}

int definicoes_obscurecimento(const int law)
{
	vector<vector<double> > par;
	
	Light_Curve teste(par);
	
	return teste.npar_obs_borda(law);
}

void show_par(ostream& out, vector<vector<double> > &par)
{
	vector<vector<string> > par_def (4);
	
	par_def[0].push_back("i = ");
	par_def[0].push_back("M1 = ");
	par_def[0].push_back("M2 = ");
	par_def[0].push_back("P = ");
	
	par_def[1].push_back("f2 = ");
	par_def[1].push_back("grav_coef = ");
	par_def[1].push_back("limb_law = ");
	par_def[1].push_back("nspot = ");
	par_def[1].push_back("T_pole = ");
	par_def[1].push_back("ffac = ");
	par_def[1].push_back("fd = ");
	
	for(int i = 0; i < par[2].size(); i++){
		ostringstream str0;
		str0 << "x" << i << " = ";
		par_def[2].push_back(str0.str());
	}
	
	par_def[3].push_back("s_theta = ");
	par_def[3].push_back("s_phi  = ");
	par_def[3].push_back("dst = ");
	par_def[3].push_back("dsp = ");
	par_def[3].push_back("Is = ");
	
	out << "# Tiago Ribeiro - LinuxC++ - 2006\n#";
	
	for(int i = 0; i < par.size(); i++){
		for(int j = 0; j < par[i].size(); j++){
			if( i < 4 ) out << par_def[i][j];
			else out << par_def[3][j];
			
			if( i == 0 && j == 0) out << par[i][j]*180/M_PI << ' ';
			else out << par[i][j] << ' ';
		}
		out << "\n#";
	}
	out << "fase  fluxo  err\n";
}

double S(vector<double> & map)
{
	double s_ = 0.0;
	double soma = 0.0;
	
	for (int k = 0; k < map.size(); k++) {
		double gx = map[k];
		s_ -= gx * gx * log(gx * gx);
		soma += ( gx * gx );
	}
	s_ /= soma;
	
	return s_;
	
	
}

double S_gauss(vector<double> & map, vector<double> & defmap)
{
	double s_ = 0.0;
	double soma = 0.0;
	double somadef = 0.0;
	for (int k = 0; k < map.size(); k++) {
		double gx = map[k];
		double gdef = defmap[k];
		
		//cout << gx * gx * log( (gx * gx) / (gdef * gdef) )  << " " << gx << " " << gdef << endl;
		
		s_ -= gx * gx * log( (gx * gx) / (gdef * gdef) );
		soma += ( gx * gx );
		somadef += ( gdef * gdef );
	}
	s_ /= soma;
	s_ += log( soma / somadef );
	
	return s_;
	
}